function [T, E, ET, TE] = build_fem_mesh(V,T)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% function [T, E, ET, TE] = build_fem_mesh(V,T)
% generate edge information for the simplest mesh V and T, secondly find 
% the neighbouring informations for simple mesh composed of V and T.
% By the way, update the vertices' index into counter clock-wise.
%
% Actually, it is only executed in the preprocessing, so the efficiency
% is not the main reason. 
%
% July 14th, 2012.
%
%  Author: Dr. Xian-Liang Hu
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% the areas of all triangles.
x21 = V(T(:,2),1) - V(T(:,1),1);
x31 = V(T(:,3),1) - V(T(:,1),1);
y21 = V(T(:,2),2) - V(T(:,1),2);
y31 = V(T(:,3),2) - V(T(:,1),2);
areas = 0.5*(x21.*y31 - x31.*y21);
% promising the index of vertices is anti-clock
pos = find(areas<0); % find the anti-direction triangles
tmp = T(pos,2); % and change their position;
T(pos,2) = T(pos,3);
T(pos,3) = tmp;
%areas(pos) = -1*areas(pos);

% first, reserve spaces for E ET TE .
nt = size(T,1);
E = zeros(nt*3,2);  posE = 0;
ET = zeros(nt*3,2); 
TE = zeros(nt,3);
%  tic;
for i = 1:nt
    for j = 1:3  % check over it's three edges
        v1 =  T(i,mod(j,3)+1);
        v2 =  T(i,mod(mod(j,3)+1,3)+1);
        edge = [min(v1,v2) , max(v1,v2)];
        edgenum = [];
        if posE > 0 % have to search the array E to check if it has existed
            edgenum = find((edge(1) == E(:,1))&(edge(2) == E(:,2)));
        end
        if isempty(edgenum) 
            posE = posE + 1; 
            E(posE,:) = edge;
            edgenum = posE;
            ET(edgenum,1) = i; % the eg is new, I'm his left triangle.
        else  % the eg has exist, then I'm the right triangle of this edge;
            ET(edgenum,2) = i; 
        end
        TE(i,j) = edgenum;
    end
end
E = E(1:posE,:); ET = ET(1:posE,:);
% bdr = find(ET(1:posE,2)==0);


end
 
%%%%%%%%%%%%%
% E = get_edge_list(T)
%
% function E = get_edge_list(T)
%     np = max(max(Elements));
%     ip1 = T(:,1); ip2 = T(:,2); ip3 = T(:,3);
%     A = sparse(ip1, ip2, -1, np, np);
%     A = A + sparse(ip2, ip3, -1, np, np);
%     A = A + sparse(ip3, ip1, -1, np, np); 
%     AA = tril(A + A');   % treat the boundary edge, which is not symmetric
%     [i,j,s] = find(AA); 
%     E = [i j];
% end

% function idx = sort_edge(V,T,E,ET,edges,k)
% 
% % find the first ring of vertex k
% nedges = size(edges,1);
% vtx = zeros(ndeges,1);
% for eg = 1:nedges
%     edge = edges(eg);  % global index of the edge
%     vtx(eg) = E(edge,1);
%     if(vtx(eg) == k)
%         vtx(eg) = E(edge,2);
%     end
% end
% 
% % check if it is a boundary vertex
% bnds = find(ET(edges,2)==0);
% if(isempty(bnds))
%     
% else
%     
% end
% 
% 
% % and then sort them in anti-clock turn
% idx = 1:nedges;  
% % idx = 1:2*nedges-1;
% % idx(nedges) = 1;
% % for k = 2:nedges  % sort the other edges according to ET
% %     
% % end
% end